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Abstract 



In this report a lumped transfer function model for High Pressure Natural Gas Pipelines is derived. Starting 
with a partial nonlinear differential equation (PDE) model a high order continuous state space (SS) linear 
model is obtained using a finite difference method. Next, from the SS representation an infinite order transfer 
function (TF) model is calculated. In the end, this TF is approximated by a compact non-rational function. 



1 



1 Introduction 



In this report we investigate the problem of the representation of a high pressure gas pipeline by a compact 
non rational transfer function model. This model is used to simulate mass flow and pressure in a small high 
pressure pipeline, and although this is a simple model with few parameters, it seems to have an accuracy 
comparable to the SIMONE'S' simulator. Since this kind of models are suitable to control design and are 
well understood by control practitioners, it is our intention to apply them to gas leakage detection and gas 
network control. 



2 STATE-SPACE DISCRETE-IN-SPACE MODEL 

The gas dynamics within the pipes is represented by a set of partial differential equations (PDE). If we 
neglect the viscous and the turbulent effects of the flow and assume small temperature changes within the gas 
and small heat exchanges with the suiToundings of the pipeline, it can be described by the one-dimensional 
hyperbolic model 

dq{i,t) ^ dpii,t) _ f,c^ q^{l,t) 

dt di 2VAp{£,t) (^^ 

dpji, t) _ dqji, t) ^' 

dt ~ ~7 d£ ' 

where I is space, t is time, p is edge pressure-drop, q is mass flow, A is the cross-sectional area, V is the 
pipe diameter, c is the isothermal speed of sound, and /c is the friction factor. 

In this research we linearised model ([T]) around the operational levels {pm.{£),qm) , where we assume a 
constant flow rate, and from the first equation of ([T|l 



pM = \lplieo)-^qU^-io). 



Hence we set p{i,t) = Pm{£) + t) and q{i,t) = qm + ^qi^^t), where Ap{i,t) and Aq{i,t) are 

deviations from the pressure/flow operational levels, respectively. Then ^ ^ ^ ^ — ^ — 



-.2 



p{£,t) p^{i) + Ap{i,t) 



+ 2^i^Aq{i, t) - -i^Ap{i, t) 



Pmii) Pm{i) ' ' Plii) 

The third term may be neglected since the distribution networks operate at very high pressure, ca. 80 bar. 
Then we substitute the remaining in the first equation 

Assumming small oscillations, Aq{£, t) ^ 2Aq{£, t), we may have {qm + 2Aq{£, t)) q{i, t) and obtain 
the following linearized model: 



dt di 
dp{i, t) _ c2 dq{£, t 

dt ~ ~A di 



.2 .^ (2) 



2 



where 

Next, decompose the pipeUne into sections £j = i = 1, 2, . . . , iV, where £q = 0, ij\r = L and L 

is the length of the pipeline. We assume the massflow to be the same in each section and accordingly define 
the following notation: 

qo{t) = q{0,t) 

qi{t) = qii,t), l,_i<i<ii, z = l,2,...,iV 
qN+l{t) = q{L,t) 

Pi{t) = p{ii,t), i = 0,l,...,N. 

Making 

iii,t) - iii.ut) 



d-ii,t) 



de 



i = l,2,...,N, (5) 



we can now approximate the hnearized PDE Q by the following discrete-in-space model 
A 

= [Pi~-i(t) - Pi{t)] - 2aqi{t), 



c2 



i = l,...,N 

j = l,...,N + l. 



(6) 



where 



Ai = i,+i-ii = ^, i = 0,...,N-l. (7) 
The pipe can then be described by the following state-space model: 

A A 

XN+l+jit) = -^fjif) - ~ '^OiXN+l+j{t) 

yi{t) = xi{t) 

y2{t) = XN+l{t), 

where i = 1, . . . j = 1, . . . N and also 

u{t) = [ qo{t) qN+i{t) f = [ ui{t) U2{t) f 

x{t)=[po{t) ••• pN{t) I qi{t) ■■■ qN{t)f (9) 
y{t) = [ Po{t) PN{t) f = [ yi{t) y2{t) f. 
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In matrix notation: 



xit) -- 

y{t) - 

Partition A as: 
A 

where 



Ax{t) + Bu{t) 
Cxit). 



^11 


Ax2 ' 


^21 


A22 _ 



A 



12 



A21 

A22 
B 







{Ar+l)x(Ar+l) 

n 

AM 



c 

'am 



A_ _A_ 




-2alN 
^2 



AAi A^e 



AAe 





A_ _'a_ 

Ai ~A£ 



€ ]R(7V + 1) X N 



AAi 



[ ei I -CAT+i 



C = [ ei I CN+i 



(10) 



(11) 



(12) 



(13) 



(14) 

(15) 
(16) 
(17) 



where e, is the i*'* vector of the canonical orthonormal basis, i.e., a vector with the the i*'^ component equal 
to one and the others equal to zero. 



3 Spectral analysis of A 

In order to leain more about the system ([3]), we analyse the spectrum of matrix A. Therefore the following 
theorem: 

Theorem 1 The eigenvalues of A defined in (I12I )- (I15I) are 

Ao = (18) 



^ = - ■'^-7'" ±n/ I 2-- sin I , , 1 1 - ( ■'^"7" 1 , k = l,...,N (19) 
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Proof: In the Appendix lAl it was proven that matrix A is equivalent to ^ up to a similarity transformation. 
Consequently they have the same eigenvalues and using (1124b . we have: 



And this is equivalent to s = and det 
value, that is, Aq = 0. 



A) = s det 


sIn - All 


-Ml 


sIn - All 


-A12 


-A21 


sIn - A22 



-A 



12 





sIn - A22 

0. Therefore, A has a zero eigen- 



From Fact 2.13.10 in |[T] pp. 62-63], we have that for arbitrary matrices A, B, C and D G H^^^ such 
that AB = BA then 



det 



A 


B 


C 


D 



det(DA - CB) 



Thus, if we take 



A 
B 
C 
D 



sljy - All 

-A12 

-A21 

sIn - A22 



we see that (s/jv — ^ii)(— ^12) = {—Ai2){sIn — An) AB = BA because slj^ — An is a diagonal 
matrix. Consequently 



det 



sIn - All 


-A12 


-A21 


sIn - A22 



det{{slN - A22){sIn - All) - ^21^12) 



Given that An = OnxN and A22 = —2alN, then 



det 



sIn - An 


-A12 


-A21 


sIn - A22 



det ((s^ + 2as)lN - ^21^12) 
det ((s^ + 2as + a^)lN - ^21^12 - o'^In) = det ((s + a)'^lN - A21A12 - o^In) ■ 



If we define the following the change of variable: 

5 = (s + a)2 
then we can write 



det 



sIn - An 


-A12 


-A21 


sIn - A22 



det {SIn - (^21^12 + o^In)) 



(20) 



(21) 



From this equation, the eigenvalues of A2iAi2+a'^lN that we denote by A (^21^12 + cP'In) , are the values 
of 5 = (s + a)^ that also set det(s/Ar -^) to zero. From (EJ, A (^21^12 + a'^Iw) = (MA) + a)^, where 
A(j4) denotes the non zero eigenvalues of A. 



From the eigenvalues properties, 

A (^21^12 + o'^In) = A (^21^12) + a^- 



(22) 
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The product A21A12 is 



^21^12 



(-Y 



-2 


1 


•• 


• 





1 


-2 


1 ■ ■ 


■ 








1 


-2 ■• 


. 











■• 


. -2 


1 








•• 


• 1 


-2 



Using Fact 5.10.25 in HI pp. 200] 
A(^i^2) = -2(^)'(l 



cos 



Then 



N + 1 



1 — cos 



, k = l,...,N 



kir 



N + l 



+ a , k = 1, , 



are the values of 5 = (s + a) that set the characteristic equation of A to zero. 
Consequently 

iAiA) + af = A{A2iAi2 + 



and from (l20l)- (|24l) the eigenvalues of A are A{A) = -a± J -2 ( -^}\ I 1 - cos 



N + l 



That is 



A(^) 



(23) 



(24) 



Recalling the definition of a in equation ([3]l, we have: 

fcC^Q ' 



and this completes the proof. 



kiT 



c 

A^'™ V2(iV + l) 



\AVAP„ 



k = l. 



The asymptotic case of the nonzero eigenvalues is reported in the following corollary 
Corollary 1 If N 00 then the eigenvalues of A are 
Ao = 



4VAP„ 



kn 



(fcC Qm 
4VAP^ 



k = l,2,... 



□ 



(25) 
(26) 
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where L is the pipe length and the time that a mass pressure takes to cross the pipeline, between its 
boundaries, at a constant speed c. 



Proof: lim Aq 
Since Ai is given by 

-4 



is trivial 



then 



sm 



kn 



2{N + 1] 



cN 



sm 



2{N + 1] 



Taking the limit when N 
and, consequently, 



lim X±k 



cN . 
oo, lim — — sm 

TV-i-oo L 



kir 



2{N + l)J 2L 



-kn 



Given that 



lim Xzrk 



AVAR, 
and this completes the proof. 




□ 



The zero eigenvalue means that there is an integrator in the pipeline model. 

It has associated an eigenvector vq, which defines a direction in the state-space where the pipeline behaves 
like a pure integrator. The following lemma gives the value of this eigenvector: 



Lemma 1 Consider A as in (fT2l) - ([T5l) . Then 

Vo = ei + 62 H h CN+l 



(27) 



where Ci is the i^^ vector of the canonical orthonormal base in IR^^"*"^, is the eigenvector associated to the 
zero eigenvalue, i. e., Avq = 0. 



Proof: Let us denote vq as 

V 

where Oat is the zero vector in and 



v = ei + e2 + --- + SN+i G K,^"^^ 
where Ei is the i^^ vector of the canonical orthonormal basis in Using (fTTI ) and 

Aiiv 



Avn 



A21V 



(28) 



(29) 



(30) 
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Given that An = 0(7v+i)x(7V+i) then Auv = 0. From ([14] ) and 



A21V 



' A 


A ' 


M 


Ai 


A 


A 


A£ 


Ai 


A 


A 




Ae . 



On 



and this completes the proof. 



□ 



Remark 1 In Lemma |7] to prove the existence of the eigenvalue associated to the zero eigenvalue we only 
used the submatrices An and A12. Therefore, we can say that this eigenvalue is generated by these sub- 
matrices. The nonlinearity of the model is only expressed by matrix A22- For this reason, if we decompose 
the full nonlinear model into a linear subsystem in cascade with a nonlinear one, the zero eigenvalue would 
appear in the linear subsystem indicating the presence of an integrator in the full model. Also, An and A12 
depend neither on pm nor on qm- 



4 Transfer functions characterisation 



We determine the transfer function. Recall that the massflow at the boundaries were chosen to be our inputs 
and the pressure at the boundaries our outputs. 



To start, we apply the Laplace transform to (fTOl ) and obtain: 

Y{s) = C{sI-A)-^BU{s) 



rj-i rj-i 

where y(s) = [ Yi{s) Y2{s) ] and = [ Ui{s) U2{s) ] , and denotes the Laplace trans- 



form of f{t). That is 



Also: 







. Us) _ 
















. Y2{s) _ 





Gn{s) Gi2{s) 

G2l(s) G22{s) 



C2 



U2{S) 



{sI-A)-^ [ B 



Bo 



Ui{s) 

U2{S) 



Ci{sI-Ay^Bi Ci{sI-A)-^B2 
C2{sI-A)-^Bi C2{sI-A)-^B2 



Ui{s) 
U2{s) 



(31) 
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and also 



Gii(s) 

G22{S) 
Gl2(s) 
G2l(s) 



Y2{s: 



U2{s)=0 



U2{S) 

Yi{s) 



Ui{s)=0 



U2{S) 

Y2{s) 



Ui{s) 



Qo{s) 
Pn{s) 



Qn+i{s) 
Pois) 



Ui{s)=o Qn+i{s) 
Pn{s. 



Qn+i{s)=0 



Q()(s)=0 



(32) 



U2{s)=0 



Qjv+i(s)=0 



4.1 Transfer function G 
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When we select this transfer function, it means that we are interested in the transfer function between the 



pressure and massflow at the intake node, i.e 



Guis) - 

and hence: 

Guis) = CiisI-Ar^Bu 



Ui{s) 
Pois) 



Qois) 



when U2is) = that is: 
^i(^) 



i(s)=0 



Ulis) 



(72(s)=0 



(33) 



where Bi and Ci are the first column and first row of B and C in ([T6l)-(fT7l). respectively. Guis) is a 
rational function whose poles are the eigenvalues of A. 



The following theorem states the zeros of this transfer function. 
Theorem 2 The zeros of Gu is) are 



Z±k 



AVAPrr 



2— sin 



(2fc-l)7r\V ffcC^Qn 



Ai \2i2N + 1)J J \4VAP, 



, k = l,...,N (34) 



Proof: To do this we recall the result from ||2j pp. 284], we have that the zeros of the transfer function ([33 
are the zeros of the following polynomial 



(35) 



SI{2N+1) - A 


-Bi 


Ci 
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RecalUng that Ci = ei G R^^iv+i) ^nd Bi = ;;^ei]R(2^+^) (see equations ^-^), we have: 



^ 



SI{2N+1) - A 


-ei 


^1 






(36) 









-1 






sIn+1 - All 


-A12 











-A21 


sIn - A22 









1 









^ 



(37) 



using the definition of An and A22 in (fT2l ) and ([13] ) 









-1 








-^12 











-^21 


(s + 2a)lAr 









1 • • 









0. 



(38) 



We develop this determinant first along the last column and next along the last row, and obtain: 



sIn 


-A12 


-A21 


(s + 2a)lN 



(39) 



where A12 denotes matrix ^412 without the 1st row and A21 denotes matrix A21 without the 1st column. 
Next, as sIn [—Ai2^ = (—^12) sIn, we apply again Fact 2.13.10 in lH pp. 62-63], which states that for 
the arbitrary matrices A, B, C and D € IR^^^ such that AB = BA then 



det 



A 


B 


C 


D 



det(DA - CB) 



and 

-L I {s + 2a) = k(- + 2"K^-^2i^2| 

= \{s^ + 2as + a^)lN - {A21A12 + 

= {s + aflN-{A2iAi2 + a^) 
We do the same change of variable as before 

S = {s + af 

and then can write: 



a 



sIn 


-A12 


-A21 


(s + 2a)lN 



\SIn - {A21A12 + 



(40) 

(41) 
(42) 

(43) 
(44) 



Next, we calculate the spectrum of matrix (^21^12 + oP') , that is: 



A (^21^12 + a^) = A (^21^12) + 



Now, we calculate the product: 



^21^12 



A 

M 






' M 

















'At 



AAi 







AA£ 
AAI 








Then 



■At 



-1 


1 


•• 


• 





1 


-2 


1 •• 


• 








1 


-2 ■• 


. 











■• 


. -2 


1 








•• 


• 1 


c 



N 



A(A2iii2) = (^)'A(M^) 



where 



-1 


1 


•• 


• 





1 


-2 


1 •• 


• 








1 


-2 ■• 


. 











■• 


. -2 


1 








•• 


• 1 


r 



Ny.N 



From 141 pp. 72] 



A (Mat) = -2 + 2 cos ( ^^^^^^ 1, A: = l,2,3,, 



,iV. 



Having that 

A (^21^12 + O?) 



-2(^)'|l-cos 



^■'^oi^y (l-cos 



\Al) 



{2k - l)7r 
2Af + 1 
{2k -1)tt 
2N + 1 



+ 
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then from (|43l) 



z±k 



1 — cos 



{2k - l)7r 
2N + 1 



/cC Qm 

~4VAPm 

fcC Qm 

and this completes the proof. 



\Ae) ^™ V2(2iV + l)y \AVAPr 



im 
'~'m 



{2k - 1)tt 
M ™ V2(2iV + 1) 



c 
sm 



\AVAP, 



2n ^ 2 

m 



k = l. 



□ 



The following corollary resolves the asymptotic case. 
Corollary 1 If N ^ oo then the zero of Gu{s) are 



z±k 



4VAP„ 



2Th 



\4:VAP„ 



k = 1,2, 



(49) 



where L is the pipe length and the time that a particle of gas takes to cross the pipeline between its 
boundaries, at a constant speed c. 

Proof: Since is given by 

-4 



then 



■ sm 



{2k - 1)tt \ cN f {2k - l)-n 
— sin 



M \2{2N + l)) L V2(2iV + l) 
Taking the limit when N ^ oo, 



cN 
lim — — sm 



{2k - 1)tt 



cN {2k - 1) TT , 
lim — — — — — lim 



sm 



{2k - 1)tt 
2(2iV + 1) J c(2/c-l)7r 



N-^oo L \2{2N + l)J N~>oo 2{2N + 1)L N^oo (2/c - 1) vr 
and, consequently. 



4L 



2{2N + 1) 



lim z±k 



'4VAPm 



±J1 



c(2A;-l)7rV ( fcC^Qr 



2L 



\4:VAP„ 



k = h2. 



L 



Given that Ta = — then 

c 



lim z±k 

iV-5-oo 



fcC^Qn 



AVAPn 

and this completes the proof. 



±J1 



(2fc-i)^ y f fcc^Q, 



2Td 



\AVAP„ 



k = 1,2,... 



□ 
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Corollary 3 Transfer Junction Gu has the following form: 



KGf[(^ + s(- + — )+l 

^ZkZ^k \Zk Z^k ' 



oo 

S 

k=l 



where Zk = a + j/Sk and z^k = a — jl3k, as well as Xk = a + jbk and X-k = a — jbk, as defined in 
Corollary\2\and Aj are defined in Corollary\J\ 

Proof: From Corollary [2] and Corollary [H we can write: 

On = ^ Z ^ <5 1) 

s \ ^ / s 



s 

k=l 



and expression (1501 ) follows immediately, after calculating the products: 
' + 1 V + 1^ and f ^ + 1 1 ( ^— + 1) . 



Zk J \—z-k ) \—Xk J \—X_k 



To complete the transfer function characterisation we need to compute the gain Kq. 
Theorem 3 Consider Bi, tite first column of B defined in (I16I ). written in the base 



□ 



{V-N, ■ ■ ■ ,V-i,VQ,Vi, . . .,Vn} 

where Vi, i = —N, —1,0,1, N, are the eigenvectores of A. If 'do is the component of Bi along vq 
then 

Kg = ^0- 



Proof: Denote the zero eigenvalue of A as Xq and Xi, i = —N, ■ ■ ■ — 1,1, . . . , N the complex eigenvalues, 
i.e. A_j = A* where •* means the conjugate eigenvalue. Since all eigenvalues have multiplicity one there are 
2N+1 independent eigenvectors Vi, i = —N, . . . , N, respectively associated to each eigenvalue Aj. Thus Vi 
and (s — Xi), i = —N, —1,0,1, N are, respectively, the eigenvectors and the eigenvalues of si — A. 
Consequently, 

{SI-A)-'v^ = -^Vi (52) 
s - Xi 

i.e., — , i = —N, . . . , N are eigenvalues of (sl — A)^^ associated to the eigenvectors Vi. Then we 

s - Xi 
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can write: 



A 



1 



s - X-N 




S — A-AT+i 



s - Ac 







s - Xn 

Also define the similaiity matrix T, considering the 2N + 1 independent eigenvectors Vi'. 

T = [ V-N ■■■ V-i Vo Vi ■■■ VN ] 

And we can write: 

{SI-Ay^ = TAT-\ 
If we decompose Bi into directions Vi, i. e., 



Bi = -d-NV-N H h '&-1V-1 + -i^o^^o + H 1- "^nvn = T 

then we can express the transfer function Gii(s) as 











Gn(s) = Ci{sl - A)-^ Bi = CiTAT-^T 




= CiTA 













N 



-V-N + 



- X-N 

After multiplying (|57] ) by s one obtains: 

■&-NS 



^ V-1 H Vo H r- 

s — X-i s s — Ai 



^1 , , 
vi-\ h 



s - Xn 



VN 



Kg = Cilim. 

s->0 \S — A_7V 

= '&0C1V0 = 'do, 



-d-is '&0S ^ -dis T^NS 

V-N -\ \ ^ — v-i H Vo H ^vi -\ \ — 

s — X-i s s — Ai s — Xn 



(53) 



(54) 



(55) 



(56) 



(57) 



VN 



□ 



because Ci = ej and vo = ei+e2 + - ■ ■ + eN+i- If we knew all the eigenvectors we could straightforwardly 
determine i9o- But, only vo is known and it is not so immediate to compute ??o- The next lemma is of good 
help to solve this problem. 

Lemma 2 : If A ^ IR"^" is a singular matrix, vq its eigenvector associated to the zero eigenvalue and 
A^vo = 0, then vq is orthogonal to the remaining eigenvectors vi, i 7^ 0, of A. 

Proof: Let Vi with i ^ Q the eigenvector of A associated to the eigenvalue Aj ^ 0. By the eigenvector 
definition 



Avj = XjVj =^ Vi = —Av. 



(58) 
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Now, using this equation, we compute the internal product between Vi and vq, 

T 

..T„ 



{vi,vo) = vl Vq = ( ^^■W) 



1 



A, 



vfA^vo = 



from which conclude that vq and Vi are orthogonal. 



Corollary 4 Consider A as defined in (fT2l) -([T5]). Its eigenvectors Vi, i = —N,...,—l,l,. 
orthogonal to vq. 



(59) 

□ 

. , A^, are 



Proof: Recalling the definitions of vq and An in equations (1281) and (1121) . respectively, then: 



Computing 



Aj^v 



4^ 
^11 


4^ " 

^21 


. ^12 


^22 . 



c 







^0 







(Ar+l)x(7V+l) 



>1A£ 
^A£ 



AM 




^A^ 
AM 



T 



AM 




^A£ 

AM 




^A£ 
AM 



AM 

r3. 



C 



AM AM 





' 1 " 




1 




1 




1 




_ 1 _ 



AM AM 



AM AM 



Then by Lemma |2] we have the expected result. 



c2 

Corollary 5 Kq = —ry ■ 
Proof: Decompose Bi as 

where is the orthogonal projection into w operator and denotes the orthogonal complement of 
From the orthogonality condition between vq and Vi,i = —N,...,—l,l,... ,N, 

PvqBi = doVQ. 

Given that 

PvqBi = {vIvqY^vIBiVq, 

then 



^0 = {vEvo)-\'Eb 



A{N + l)M' 



Replacing by ^ we find 

° A{N + l)L- 

N 

When iV oo ^ — > 1, 

iV + 1 



Corollary 6 Transfer function Gu has the following form, according to Corollary\3\ 



/I 1 



+ s\ — + +1 



where Zk are defined in Corollary^and Xk are defined in Corollary\l} 



4.2 Transfer function G 
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Next, we determine the transfer function 

Y2{s) PNis) 



G 
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U2{s) Qn{s) 
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when Ui{s) = 0. That is the ratio between the pressure and the massflow at the offtake node. From (|3T]) : 

G22{s) = C2{sI-Ar^B2, (68) 
where C2 is the second column of C and B2 is the second row of B. 

Similarly to what happens with Gu{s), ^22(5) is a rational function whose poles are the eigenvalues of A. 
Following the same methodology as for Gii(s), we would like to calculate the zeros of G22(s) in order to 
investigate pole-zero cancelations. 



Theorem 4 The zeros of G22{s) are 

fcC^Q ^ 



Zk 



AVAR 



2N + 1 



\AVAP„ 



k = l. 



(69) 



Proof: The proof is very similar to the one of Theorem|2l Again, we recall the result from E pp. 284] that 
states that the zeros of the transfer function (l67l ) are the zeros of the following polynomial: 



(70) 



SI{2N+1) - A 


-B2 


C2 






Recall the definitions of C2 = cn+i G M^^^^^^ and B2 
exactly as for Theorem |2] 



GN+l 



and the proof follows 



We have: 



SI(2N+1) - A 










<^ 



sIn+1 - All 


-A12 


-A21 


sIn - A22 







1 







_0_ 

T 









From the definition of An and A22 defined in (fT2l) and ( fT5l ) 

= ^ 


















sIn+1 


-A12 













_ -A21 


{s + 2a)lN 






= 0. 















••• 


0\1 ••• 








(71) 



(72) 



(73) 



(74) 
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We develop this determinant first along the last column and next along the last row, and obtain: 



sIn 


-A12 


-A21 


(s + 2a)lN 



(75) 



where A12 denotes matrix A12 without the 1st column and A21 denotes matrix A21 without the 1st row. 
Next, as sIn (—^12) = (—^12) sIn, we apply again Fact 2.13.10 in HI pp. 62-63], that states that for the 
arbitrary matrices A, B, C and D € JR^^^ such that AB = BA then 



det 



A 


B 


C 


D 



det(DA - CB) 



and 



sIn 


-A12 


-A21 


{s + 2a)lN 



s{s + 2a)lN - A2iAi2\ 

(s^ + 2as + a^) In - (^21^12 + a^) | 

{s + af In - {A2iAi2 + a'^) . 

We do the usual change of variable 

S = {s + af 

and then can write: 

= |5/7V- (^i^2 + a^)| 
Next, we calculate the spectrum of matrix [A21A12 + a^) , that is: 

A (^21^2 + a^) = A (^1^12) + 
Now, we calculate the product: 



AoaA 



21^12 



_A_ 

Ai 
A_ 

Ai 









A_ 

' A£ 














A_ 

Ai 




_A_ 
Ai 



AAi 




AA£ 
AAi 



AAi 








-_C_\2 

-Ai' 



-2 


1 


•• 


• 





1 


-2 


1 •• 


• 








1 


-2 ■• 


. 











■• 


. -2 


1 








•• 


• 1 





€ IR 



N 



(76) 

(77) 
(78) 

(79) 
(80) 



A^i 



AAi 





c2 



AAi J 
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Then 



I SIJV -^21^2 I = (^)2A(M7v) 



where 



Mm 



-2 


1 


•• 


• 





1 


-2 


1 •• 


• 








1 


-2 ■• 


. 











■• 


. -2 


1 








•• 


• 1 


-1 



€ IR 



NxN 



(81) 



From H pp. 72] 

A(Mjv) = 



-2 + 2cos( ^^^^ ) , A: = 1,2,3,. ..,iV. 



(82) 



and from here the proof follows exactly as for the calculus of the zeros of Gii(s), and we can see that are 
the same. 



Likewise follows for the asymptotic case. As we can see from the definition of the transfer function (1331) 
and (|68] ) as well as from the definition of the Cj, i = 1, 2 we have 

Gii{s) = -G22{s) (83) 

Therefore, its zeros will be necessarily coincident. 



□ 



4.3 Transfer function G12 

Acccording to (|3TI ). consider now the transfer function 

_ Y,{s) _ Pojs) 

"""^'^ -TW)- Q^) ^''^ 

with Ui{s) = Qo{s) = 0, or equivalently: 

Gi2{s) = Ci{sI-A)'^B2, (85) 
where Ci is the first column of C and B2 is the second row of B. 



Theorem 5 The transfer function Gi2{s) has no zeros. 

Proof: According to IH pp. 284], we have that the zeros of the transfer function (l84l ) are the zeros of the 
following polynomial: 



SI{2N+1) - A 


-B2 


Ci 
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and recalling that Ci = d e IR^^TV+i) = -eAr+i]R(2^+^\ we have: 



sI{2N+l) - A 










= ^ 
















sIn+1 - All 


-A12 






1 









-A21 


sIn - A22 






1 









And from de definition of An and A22 in (fT2l ) and ([T5] ) 

= 4^ 
















sIn+1 


-A12 




1 









-A21 


{s - a)lN _ 




1 ••• 









We develop this determinant first along the last column and next along the last row, and obtain: 









•■ 


• 








7 








• 








s 





•■ 


• 








-7 


7 














s 


•• 


• 











-7 


7 














■■ 


s 

















• -7 


7 











•• 


■ 


s 














• 


-7 


7 


Q 





•• 


• 








(s + 2a) 








• 








-Q 


Q 


•■ 











(s + 2a) 














-Q 


Q ■ ■ 














(s + 2a) • • 














• -Q 


Q 











• 


(s + 2a) 











•• 


■ 




Q 











• 





(s + 2a) 



= 



(88) 
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We don't worry about the signs of the cofactor, since our aim is to determine the zeros of the determinant. 
Now, we develop this determinant first along the first line, and we obtain: 



7 



s 





• 


• 








7 





• 











s 


• 








-7 


7 














s 












• -7 


7 











• 


■ 


s 











■ 


-7 


7 


Q 





■ 


• 














• 








-Q 


Q 


• 









(s + 2a) 














-Q 


Q ■ 












(s + 2a) • ■ 














• -Q 


9 











(s + 2a) 











• 


■ 


-Q 


Q 








• 





(s + 2a) 



Next, we develop along column- A^: 



s 





•• 


• 





7 





■ 











s 


•• 






-7 


7 














s 









• -7 


7 











•• 


• 


s 








■ 


-7 


7 


Q 





•• 


• 











■ 








-Q 


Q 


•• 






(s + 2a) 














-Q 


Q ■ ■ 


• 








(s + 2a) • • 














•• 


• -Q 


9 








• 


s + 2a) 






Here the matrix is of dimension 2{N — 1) x 2{N — 1). 
Swap the first row of blocks with the second one: 



Q 





•■ 


• 











• 








-Q 


Q 


•■ 






(s + 2a) 














-Q 


Q ■ ■ 









(s + 2a) • • 














•• 


• -Q 


9 








• 


s + 2a) 





s 





•• 


■ 





7 





■ 











s 


•■ 






-7 


7 














s 









■ -7 


7 











•• 


■ 


s 








■ 


-7 


7 



Swap the 



first column of blocks with the second one: 









• 













•• 


■ 





(s + 2a) 













Q 


•• 









(s + 2a) • • 











-Q 


Q ■ ■ 












• 


(s + 2a) 











•• 


■ -Q 


Q 


7 





• 








s 





•• 


■ 





-7 


7 











s 


•• 










• -7 


7 











s 











• 


-7 


7 








•• 


■ 


s 



Again develop the determinant along the first row: 



{s + 2a) 





■ 










•■ 


• 








(s + 2a) • • 










g ■ ■ 












■ 


(s + 2a) 








•■ 


• -Q 


B 


7 





■ 











•• 


■ 





-7 


7 








s 


•• 










• -7 


7 









s 











• 


-7 


7 





•• 


• 


s 



Again along column-(A^ — 1) 

(s + 2a) 
(s + 2a) 



2 2 

7 Q 











(s + 2a) 



Q 








-Q Q 



•••0 ••• 

7 • • • : : s 



Therefore, we find a pattern. To write the pattern, we define: 

In = identity matrix whose first row and column- are all zeros 

Ai2 = matrix A12 without the last row 

^i2jv = matrix of order A^ and with the same pattern as A12 

A21 = matrix A21 without the first column 

^2ijv = matrix of order A^ and with the same pattern as A21 
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with this notation, we can write the determinant ([88]) as: 













{s + 2a)lN~i 



with ? = 0, it becomes 



Define an iteration as: 

1. Develop the determinant in cofactors along the first row 

2. Develop the determinant in cofactors along column {N — i) 

3. Switch the first row of blocks with the second one 

4. Switch the first column of blocks with the second one 
Then, we obtain: 







{s + 2a)/Ar_(i+i) 




^12]v-(i+i) 


slN-{i+l) 



Iterate again and obtain: 



y+2^i+2 







slN-{i+2) 


^12]v_(i+2) 


^21]v_(i+2) 


(s + 2a)/;v_(j+2) 



Also, considering = 1 we write (1881) as: 



s 





-7 








s 


7 


1 


Q 


-Q 


(s — a) 





1 












(89) 



Similarly to what we have done for the general case, we develop the determinant first along the last column 



s 





-7 


Q 


Q 


(s — a) 


1 









and next along the last row: 






-7 


Q 


{s — a) 



and the proof that Gi2(s) has no zeros is complete. 



□ 



Gi2(s) is a rational function whose poles are all the eigenvalues of A, since this transfer functions has no 
zeros. Therefore, we can write: 
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Corollary 7 Gi2{s) is given by 

+ s{ — + - +1 



where \±k i^ given by ( |26l ). 



5 Approximated Transfer functions 

In this section we propose some approximations for the models of the transfer functions. 
From Corollaries [T]-E] and Theorem [3l Gij{s), i, j = 1, 2, are meromorphic functions given by 



Gu[s) = II Kk ; ^2 TTTl ^ — 



G,M = ^\\ , (91) 

G22(s) = -Gii(s) 
Gi2(s) = -G2i(s) 



where 

2Trf 
and 



'^o = 7^ (92) 



Natural gas is highly pressurized in transportation networks in order to expedite its flow. To ensure this, 
it must compressed periodically along the pipe. This is accomplished by compressor stations, which are 
usually placed at 60 Km to 250 Km intervals along the pipeline. As a result, the frequency ujq always 
remains much greater than a. Taking this into account as well as the requirement that its factors have a DC 
gain set to 1, we define A'^, and thence can approximate Gij{s), i,j = 1, 2 by 



r Kg fip {s + af + {2k - ifu^l 



k=l 



Gi2(s) = -G2i(s) 
G22{s) = -Gii{s) 



with 



+ Ak'^ujl 
a2 + {2k - IYujI 



Kk = :^r:^. 2 - (95) 
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If we define S = s + a we can write 

Guis) - G'n(S) - g-^ 11 K, ^ 

(96) 

k=l " 

Theorem [TO] considers auxiliary functions that lead to significant simplification in the representation of 
Gii(S) and G2i(S). However, before stating Theorem [TOl we need to prove some intermediate results: 

Theorem 6 The function 



1 - e- 



/(«) = . ^-2sT, (97) 



be expanded as 



/(.)=E7^ (98) 



k=—oo 



where Uk is the residual of f{s) at s = Afc, /. e. 



.l)k 



2T, 



Since condition ( 11291 ) holds (see appendix |B]| then the expansion ( |98] ) exits and the residuals are given by 
Proof: 

Ofc = 7T^ / f{s)ds = lim (s - \k)f{s) = lim — — 



lim ^"'^^ ~ ~ Xk)Tde-'^^ _ e^^^-i _ e^^" _ e^^^ _ {-if 
Theorem 7 Thefiinction 



□ 



(99) 



OO / s 



A;=— 00,^7^0 



n ('-^"tJ 

k=—oo 



can be written as: 
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Proof: Since g{s) is proper it can be expanded as a Laurent series 



s - Xk 

k=—oo 



where bk is the residual of g{s) at s = Xk = jk — . In order to compute this residual we rewrite g{s) as 

Td 



M . . 

n {-'"y) 



, s ,. k=-M,kj^O 

g[s) = nm 



n (-i^) 

k=-M " 



The residual is then given by 



h = 77^ / 9{s)ds = lim (s - \k)g{s) = lim — — - ^^•'"^o 



JJ^ ( jni— I 



m=—M ^ ■m=fc+l 



Given that 



TT \ / vr \ n / mvr 



-J?7T,— jm 



we can rewrite bk as 



M . X 2 
TT 



bk = lim ^ 



M-i-oo / „ \ M 



=-M ^ m=fc+l ^ 
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If now we replace define £ = k — m vje get 

A / ^ N 2 i\ / / \ 2 

n rij n (^^^ 

, ,. m=l ^ '^^ v m=l ^ ^ 

h = lim — — = lim 



u O^k) n i'"n fe n 



e=-(M-k) ^ i=i ^ '^^ f=_(M-fc) ^ ''^ 

n 



1- m=l 

iim 



M ^ \ 2 

TT 



M->oo / _ \ A/-fc ^ N 

n 



1- m=l 

lim 



Af / X 2 

TT 



n (t^) n (4) 

A^/xAZ/x Ai". 

n ^^ n Oi n 

lini ; — — — — ; = iim 



M^oo k+M / \ M-k / ^ N A/^oo A/+fc / n 

(-1)="-' n {'^) n (-')-' n {%) 



j.^ fc [M-(A;-l)][M-(fc-2)]...(M-l)M 
M^oo^ ' {M + 1){M + 2)...[M +{k-l)]{M + k) 

AfToo*^" M + /c M + (fc - 1) " ' M + 2 M + 1 

, M-f/c-l) , M-(k-2) , M-1 , M 
(-1)" lim ^ ^ lim ) f . . . hm lim 



A/^oo M + k M^oo M +{k-l) M^oo M + 2 A/^oo M + 1 



k 



i-l)"!] lim 



M -(i-l) 



A/^oo M + £ 

£=1 

Since we can always make M infinitely greater than k then 

A^^OO M + i 

and, consequently, 

6fe = (-1)'= = 2Tdak 
and we conclude that 



n i-^"' 



= 



n (^-i^f ) 

fe=— oo 
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and the proof that the expansion of function /(s) in Laurent series is possible is done in Appendix iBl 

□ 

Theorem 8 The function 

l + e- 
can expanded as 



oo 



k=—oo 

where is the residual ofv{s) at s = ■ 

'''' 2Td 

Proof: This function has poles at 

vr 

\k=j{2k-l) — , /c = -oo,...,-l,0, l,cx) 
2'J^d 

We can prove that condition ( 11291 ) holds for v{s) exactly the same way we did for f{s). So, we can expand 

v{s) as 



k=—oo 



where is the residual of v{s) at s = A^, i. e. 



Cfe = -— f v{s)ds = hm [s - Xk)v{s) = lim — — ^r;— 

2j7r Jcfe ^-^^fc 1 + e '^-'<'* 



j(2fc-l)7rrp 



□ 



Theorem 9 The function 



k=—oo 

can be written as: 



n -^-(^'-^'^ 



w{s) = 55 ; — (101) 
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Proof: 

Since w{s) is proper it can be expanded as 



oo , 



s - Xk 

fe=— oo 

TT 

where is the residual of w{s) at s = = j{2k — 1)——. In order to compute this residual we rewrite 

w{s) as 



M 

n 

k=-M+l 


(-i(2A: - 




M 

n 


s-i(2A;- 





?i;(s) = lim 



The residual (i^ is then given by 

M 



dk = — — f w{s)ds = lim {s — Xk)w{s) = lim ^^"""^ 



2j7r /c^ s^Afe J\/->oo ^1 / TT \ -j^ 



n fi 6(A^-"^)f) 



m=— A/+1 m=k+l 

Given that the numerator can be expressed as the product of two complex factors we can write dk as 



M / ^ 2 

vr \ 



dk = lim ™ ^ 



n (2--!)^ 



2Td / 



n (j{k-m)-\ H (j{k-m)-\ 

m=-M+l ^ '^^ m=k+l ^ 

If now we replace define £ = k — m we get 

M . X 2 M 



n((— )^) ^ n((— 



dfc = lim - — — = lim 



n n n 0| n 

A/ . X 2 

n p-D^ 

m=l ^ 

(-i)-'r' n (^i) n (-^ 

M . s 



lim 

A/ / N 2 



J(-I)'^' 
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and we conclude that 



OO / X 



n (^-j(2fc+i)^) 

fc = — OO 



Theorem 10 Consider the following functions 



□ 



OO ^ \ 



Ge(S) = ^" ^'^'""^ (102) 

vr 



n (s-^'=4) 



fc = — OO 



OO / X 

n -ip*-')^ 

OO ^ X 



0„(S) = 4^^^^ (103) 



fc = — OO 

FM = 1 _ ^-2T,s (104) 

^o(^) = 1 + g-2ST. (105) 



Ge(S) = 2TdF,{S) (106) 
Go(S) = 2Fo(S) (107) 

Proof: To prove the results in Theorem |7] we consider g{s) = Ge(s) and /(s) = Fe(5) and also in 
Theorem|9]we consider w{s) = Go{s) and v{s) = -Fo(S'). 

□ 

Ge(S) and Go(S) may also be written as 



OO 



„2, .2 



Ge(S) = (108) 



k=l 

OO 

k=l 

n (s^ + {2k - i)w,) 



k=0 

OO 



\2 2 



ColS) = (109) 



k=i 
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where Uo is defined in equation ( [92] ). As a result, 

°° + (2fc - ^ SGe(S) -rj 1 

OO -j^ OO -j^ 

nS2+4fc2 2 = SGe(S)n^^ 

k=l 

with Kk being defined in equation (|93] ). Now, using these equations and theorem [lO] we can write 

+ (2fc - 1)2^2 _ SGe(S) _ 

^ i^nT.^ (110) 



with 



°° n'2 -L 4^2, ,2 

^12 = 11 4^2 2 

fc=l 



(111) 



we thus can rewrite Gii(S) and G2i(S) as 
(S-a)(l-e-2ST,) 

(112) 

(S-a)(l-e-2STrf)- 

We can compute Gii(s) and G2i(s) from these equations by replacing S with s + a: 

^ ^ ^ KcT.Kn {s + a) (1 + e-2°^c.e-2.T,) 



G2i(s) 



s(l_e-2"T,g-2sTd) 

(113) 

5(l_g-2aTde-2sT,) 

Given that Gn (s) and G21 (s) were defined in such a way that 

limsGii(s) = lim sG2i(s) = lim sGii(s) 

s— >o s—^O s—>-0 (114) 

= lim sG2i(s) = Kg, 

we can rewrite ( 11131 ) as 

^ (. + a) (1 + e-2"^^e-2-^^) 

Gll(^) = ^(i_^_2.T,e-2sT,) 

(115) 

(s + a) e"'*^'* 



G2l(5) =^21 



s(l_g-2aTrfg 



2Q:Trf^-2sTrf 
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with 



K 



Kg {I 



-2aT, 



11 



0(1+6-2°^^) 



K' 



Kg (I 



-2aTd 
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a 



(116) 



6 Case study 



In this section we study a small pipeline using the lumped linear model derived above. The pipe has a 
length L = 35Km and a diameter V = 793mm. The friction factor is fc = 0.0079 and the isothermal 
speed of sound is c = 300m/sec. We considered qm = 90Kg/sec and pm = 80 bar = 8 x 10^ Pascal 
as the mass-flow and pressure nominal values. The a and Kq parameters were calculated from equations 
Q and Corollary |5l respectively, and the values of 0.0051 and 5.064 were obtained. We computed the 
frequency responses (FR) of Gii(s) and G2i(s) from equations ( |9T| ) using truncated approximations of 
order n, G'^"''(s) and G^^\s), respectively, in a frequency bandwidth, BW = [lO~^rad/sec, 1 rad/sec] . 



BV_12400_A- 



JCT_LEAK 



28.49 Km 



TERMINAL 



7.09 Km 



Figure 1 : Gas pipeline topology. 



Figures [2] and |3] display the respective Bode diagrams and compare them with their approximations Gii(s) 
and 6*21 (•5)- G'ii(s) converges very fast. The FR of a truncated approximation with two hundred poles and 
zeros had already converged to its limit in the whole frequency interval BW. Figure |2] shows that there 
ai^e no significant differences between Gii(s) and Gii(s). Consequently, Gii(s) can be substituted by its 
approximation without loss of accuracy and with a significant reduction of the computational costs. 

The convergence of G21 (s) is much slower. A two thousand order approximation didn't converge in the 
whole bandwidth BW. We can subsequently conclude that a truncated approximation of equation ( |9T| ) 
needs too many factors leading, therefore, to high order transfer functions with high computational costs. 

However, from Figures|3]and|4]one can conclude that G2i(s) can be also substituted by its approximation. In 
Figure lUthe Bode diagrams of G21 [s] and G^2i^'^^ i^) compared with a better resolution, i.e. the phase is 
restricted to the range (—180 degrees, 180 degrees]. The phase discontinuities are due to the phase-crossing 



50 
40 



10"'' 10"^ 10"^ 10"' 



rad/sec 




Figure 2: Bode plots of Gii(s) and the approximation Gii(s). 
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Figure 5: Input and output mass-flows on the both days. 



of the odd multiples of 180 degrees which are converted from -180 degrees to 180 degrees. Also notice that 
the Bode diagram of the finite approximation G^^2^'^^ i^) converges to Gi2(s) almost in the whole frequency 
BW. 

This pipeline was simulated taking two normal days operation data as the input and output mass-flows. The 
simulation was performed with the previous refen^ed SIMONE® simulator. Figure |5] shows the input and 
output mass-flows on both days 

The intake and offtake massflows are denoted by qi{t) and q2{t), respectively. 

Figures [6] and |7] compare the input and output pressures simulated with the lumped transfer function model 

Pi{s) = Guis)Qi{s) - G2i{s)Q2is) 

P2is) = G2lis)Qiis) - G22{s)Q2is) 

with the ones simulated with SIMONE® for the two days. The intake and offtake pressures are denoted 
by pi{t) and P2{t), respectively. We can see that the results are better for the first day. This is expectable, 
since the mass-flows and pressures of the second day data present stronger deviations from the nominal 
values used for the calculation of the a parameter. As a matter of fact, on the second day the quotient 
+ Q2{t)) I {pi{t) + P2{t)) has a root mean square deviation from the nominal value of about 50% 
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Figure 6: Input and Outputpressures, Pi, and P2, respectively, simulated by the lumped model (1117b (solide 
red line) and by SIMONE® using the first day data. 
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Figure 7: Input and Output pressures. Pi, and P2, respectively, simulated by the lumped model (1117b (solide 
red line) and by SIMONE® using the second day data. 
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(this deviation is of 30% in the first day). Yet, in both cases, the model has well captured the dynamics of 
the system and, therefore, it seems to be a valuable tool for gas leakage detection and gas networks controller 
design. 



A Change of variables in an integral model for a short gas pipeline 

In this section, a change in the state-space variables of the integral model is performed. The purpose is to 
obtain a simpler system matrix in order to simplify the determination of the eigenvalues of the system. 

Consider model (ITOl) with the respective matrices defined according to (fT2)) -(fT7]). 

Consider for this system the following change of variables: 

Z2 = X2+X2-\ \-XN + XN+l 

Z2 = X2 — X2 

Z3 = X2- X3 

ZN+1 = XN-XN+1 (118) 

ZN+2 = Xn+2 

ZN+3 = Xn+3 

Z2N+1 = X2N+1 

and we obtain the following realisation 



ZN^lit) = ^e^Mt)-2^^Z2N^,{t) + -^^U2it) 
U\ U\ Qml /,N 



(120) 
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ZN+i{t) 



Z2N+l(t) 



A. , ,s, fcC Qml 



ZN+l{t) 



2VAP, 



Z2N+l{t) 



ml 



(121) 



y2{t) 

2/2 (t) 



1 



+ 1 
1 



Z2{t) + 
Z2{t) - 



N 



N + 1 

and in matricial form, we have: 

z{t) = Az{t) + Bu{t) 

with 



N + 1 
1 

N + 1 



Z2{t) + 
Z2(t) - 



yit) = 


Cz{t) 








' 


OlxTV 


OlxAf 


A = 


OatxI 


All 


A12 




. OatxI 


A21 


A22 


An = 


OatxAT 







N - 1 
N + 1 
2 

ivTT 



zz{t) + • 
zzit) - ■ 



1 



N + 1 

N 
N + 1 



ZN+l{t) 
ZN+l{t) 



(122) 



(123) 



(124) 



(125) 



A 



12 



AM 

^2 



AM 




AM 

^2 



AM 

^2 



AM 



AM 

r2 



AM 



AM 




AM 
AM 



AM 
AM 



GM^^^ (126) 
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Ml 

B 
C 



A 

fcC^Q 



In, In — identity matrix N x N 
In 



2VAPr, 



[ e2 I -62 + CN+i ] G ]R(2^+^) 



x2 



1 



N 



N -1 



1 



N + l N + 1 N + 1 
1 1 2 



N + 1 
N 



jj^2x(2Ar+l) 



(127) 



Matrices A and A have the same spectrum, however its calculation seems to be easier if we use matrix A. 



B Rational expansion of meromorphic functions 

Let /(s) be a function meromorphic in the finite complex plane with poles at Ai, A2, . . ., and let (Fi, r2 , . . .) 
be a sequence of simple closed curves such that: 



• The origin lies inside each curve Ffc. 

• No curve passes through a pole of /. 

• Ffe lies inside T^+i for all k 

• lim d{Tk) = 00, where d{Tk) gives the distance from the curve to the origin 
Suppose also that there exists an integer p such that 

fis) 



lim 

k—^oo 



\ds\ < 00. 



Denoting the principal part of the Laurent series of / about the point Afc as PP [f{s); s = Xk], we have, if 
p < 0. 



fiz) = J2^Fif{z);z = X, 



k=0 



Theorem 11 Consider function f{s) 
f{s) 



. There exists an integer p such that 



lim 



\ds\ < 00. 



(128) 



Proof: This function has poles at 

\k= 3'2k^ = jk^, k 
^Id J-d 



-00, 



-1, 0, 1, 00 



vr vr 

The contours will be squares vertices at +{2k — 1)^^ ^ — 1)^^' k > 1, traversed counterclock- 
wise, which are easily seen to satisfy the necessary conditions. To see what are the terms of the Laurent 
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Figure 8: Contours 



series expansion of /(s) we need do see for which p the condition 



Re 



lim 



fis) 



sP 



\ds\ < oo 



(129) 



holds. We can partition this integral as 



sP 



\ds\ 



where 
Notice that 



^1 



s:s = {2k-l)^^+ju;, ut, a; G [-(2/c - 1)^, (2fc - 1)]-^^ 
s:s = a + j2{k-l)^^, a i, a G [-(2A. - 1)^, (2A. - 
s:s = -{2k-l)^+jio, [-(2A:-l)T5^,(2fc-l)]-^ 
s:s = a-j{2k-l)^^, ^ |, ^ e [_(2fc _ 1)-^, (2fc - 1)]^ | 



For s G Ffci, |(is| = du. 
For s G rfc2, \ds\ = —da. 
For s G Tfca, |(is| = — dw. 
For s G rfc4, |(is| = da. 



39 



Now we can write 



lim ® 



SP 



\ds\ 



+ 



-(2fe-l)^ 

(2fc-l)^ 
(2fc-l)^ 

■(2fc-l)^ 
(2fc-l)5i- 



(2fc-l)^ 



(2fc-l)5f- 



(2^-1)5% 



vr 



/((2fc-l)^+j^ 



(2fc-l)#T 



vr 



/(-(2A:-1)— 



r{2k- 






duj + I 






J~{2k- 







fU-j{2k-l) 



vr 



da 



vr 



/( (2A:-l)^+ia; 



{2fc-l)^ 
-{2fe-l)^ 



/ a + i(2fc-l) 



vr 



f{-{2k-l)—+ju 



A2k- 






duj + 






J-(2k- 







TT 

2T^ 

) — 

^^2Trf 



da 



Next we analyze the four terms of this integral 
First term 



(2fc-l)5 
■(2fc-l); 



TT 



/((2A:-1)— +ia; 

(2k~l)^ 



Given that 



dbJ 



{2fc-l)2 
-(2fc-l), 



-(2fc-l)f g-jc^Td 



(2A; - 1)^ + jo;) (l - e-(2fc-iKe-i^^rf) 

(2fe-l)f 



{2k - 1)2^ + jt^j (1 - e-{2fc-i)^e-i'^^d) 



then 



lim 

k^oo 



M2k~ 












'-{2k- 


-1) W7 





lim e-^^'^-i)! = 



vr 



duj = lim 



(2k-l)^ 
-(2fe-l)w- 



-(2fc-l)f 



= 0. 



Second term 















/-(2fe- 







(2A;-l)^ 

-(2fc-l)2T- 



g-aTdg-j(2fc-l)^ 



da 



-(2fe-l), 



(a + j(2fc-l)2f-) (l-e-2-r<^(-l)'=) 
Notice that 



da 



(2fc-l): 



(a + ji2k-l)^yil~e-^^^''{-m 



da 



-(2fc-l)2f- 



<7 + j(2fc-l)2|- (l-e-2-T.(_i)fc) 



da < 



-(2fc-l)2f- 



a + :ii2k-l)^ (l_e-2.T.) 



dcr 



On the other hand, for a < 



1 -e" 



1 < e 
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As a result 



f 



-o-Td 



da < 



f 



da 



f 



r,crTd 



Making k ^ oo, 



.0 
lim / 



a + ji2k-l)^ 



da = Mi.i 



da 



CT+j{2k-l)^^ 



Notice also that 

(2fe-l)2f- 



For cr > 



1 - e"'"™ < 1. 



da < 



(2fe-l); 



Consequently 

'•(2fc-l)#- 



/ 



Making k ^ oo again 



da < 



(T + j(2A:-l)2fj 



lim 



(2fc-l); 



a + i(2/c-l)2fj 



da = Ml 



and we conclude that 



lim 





'2Td 










/-(2fc- 







f{a + j{2k-l] 



vr 

21^ 



dcr < 2Mi < oo. 



Third term This term is similar to the first one and using the same arguments we can prove that 



lim 





'■2Td 










l-{2k- 







TT 
'21^ 



duj = 0. 



In similar way that we did for the second term we can prove that 



(2fc-l)2 
-(2/0-1), 



TT 



da < 2Mi 
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We can now conclude that condition ( 11291 ) holds for any p < 0. 

Since there exists an integer p such that 

fis) 



□ 



lim 



sP 



\ds\ < oo. 



Then, denoting the principal part of the Laurent series of / about the point A^. as PP [/(s); s = A^], we 
have, if p < 0. 

oo 

/(z) = ^PP(/(z);z = Afc). 
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